# Clears work environment
rm(list = ls())

# Sets working directory
setwd('C:/Users/Jason/Box Sync/Home Folder jdt34/733 - Maximum Likelihood Estimation/Lai & Slater 2006 replication')

# Loads data and libraries
load(file='cleanedLSdata.RData')
loadPkg(packs)

# Constructs scenarios
scenRegimesLag0 = with(lsm, cbind(1,
                                  c(0, 1, 0, 0, 0), 
                                  c(0, 0, 1, 0, 0),
                                  c(0, 0, 0, 1, 0),
                                  c(0, 0, 0, 0, 1),
                                  median(total), 
                                  mean(cap), 
                                  mean(open), 
                                  median(totalally), 
                                  0))
scenRegimesLag1 = with(lsm, cbind(1,
                                  c(0, 1, 0, 0, 0), 
                                  c(0, 0, 1, 0, 0),
                                  c(0, 0, 0, 1, 0),
                                  c(0, 0, 0, 0, 1),
                                  median(total), 
                                  mean(cap), 
                                  mean(open), 
                                  median(totalally), 
                                  1))

# Prepares colors
panelcols = hcl(h=seq(15, 375-360/5, length=5)%%360, c=100, l=65)
regcol = combn(hcl(h=seq(15, 375-360/5, length=5)%%360, c=100, l=65), 2)
regcom = combn(c('Democracy', IVs[1:4]), 2)
colcom = NULL
for(i in 1:10) {
  coltemp = as.vector(regcol[,i])
  names(coltemp) = regcom[,i]
  colcom = c(colcom, coltemp)
}

# Creates diagonal panels with descriptive statistics
# Constructs diagonal panels
panel.dem = ggplot(data=data.frame(cbind(c(0, 1), c(0, 1))), aes(x=X1, y=X2)) + 
  geom_blank() + 
  theme(panel.background=element_rect(fill='white'),
        panel.grid.minor.x=element_blank(), 
        panel.grid.minor.y=element_blank(),
        panel.grid.major.x=element_blank(), 
        panel.grid.major.y=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.ticks.x=element_blank(),
        axis.ticks.y=element_blank()) +
  annotate('rect', fill=panelcols[1], xmin=0.07, xmax=1, ymin=0.12, ymax=1, alpha=.5) +
  annotate('text', x=.55, y=.85, label='Democracy') +
  annotate('text', x=.55, y=.6, label='N = 1565 (32%)', size=2.5) +
  annotate('text', x=.55, y=.3, label='MID country-years:\n140 (22%)', size=2.5)

panel.mac = ggplot(data=data.frame(cbind(c(0, 1), c(0, 1))), aes(x=X1, y=X2)) + 
  geom_blank() + 
  theme(panel.background=element_rect(fill='white'),
        panel.grid.minor.x=element_blank(), 
        panel.grid.minor.y=element_blank(),
        panel.grid.major.x=element_blank(), 
        panel.grid.major.y=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.ticks.x=element_blank(),
        axis.ticks.y=element_blank()) +
  annotate('rect', fill=panelcols[2], xmin=0.07, xmax=1, ymin=0.12, ymax=1, alpha=.5) +
  annotate('text', x=.55, y=.85, label='Machine') +
  annotate('text', x=.55, y=.6, label='N = 1438 (29%)', size=2.5) +
  annotate('text', x=.55, y=.3, label='MID country-years:\n183 (29%)', size=2.5)

panel.jun = ggplot(data=data.frame(cbind(c(0, 1), c(0, 1))), aes(x=X1, y=X2)) + 
  geom_blank() + 
  theme(panel.background=element_rect(fill='white'),
        panel.grid.minor.x=element_blank(), 
        panel.grid.minor.y=element_blank(),
        panel.grid.major.x=element_blank(), 
        panel.grid.major.y=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.ticks.x=element_blank(),
        axis.ticks.y=element_blank()) +
  annotate('rect', fill=panelcols[3], xmin=0.07, xmax=1, ymin=0.12, ymax=1, alpha=.5) +
  annotate('text', x=.55, y=.85, label='Junta') +
  annotate('text', x=.55, y=.6, label='N = 150 (3%)', size=2.5) +
  annotate('text', x=.55, y=.3, label='MID country-years:\n28 (4%)', size=2.5)

panel.bos = ggplot(data=data.frame(cbind(c(0, 1), c(0, 1))), aes(x=X1, y=X2)) + 
  geom_blank() + 
  theme(panel.background=element_rect(fill='white'),
        panel.grid.minor.x=element_blank(), 
        panel.grid.minor.y=element_blank(),
        panel.grid.major.x=element_blank(), 
        panel.grid.major.y=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.ticks.x=element_blank(),
        axis.ticks.y=element_blank()) +
  annotate('rect', fill=panelcols[4], xmin=0.07, xmax=1, ymin=0.12, ymax=1, alpha=.5) +
  annotate('text', x=.55, y=.85, label='Boss') +
  annotate('text', x=.55, y=.6, label='N = 1172 (24%)', size=2.5) +
  annotate('text', x=.55, y=.3, label='MID country-years:\n161 (25%)', size=2.5)

panel.str = ggplot(data=data.frame(cbind(c(0, 1), c(0, 1))), aes(x=X1, y=X2)) + 
  geom_blank() + 
  theme(panel.background=element_rect(fill='white'),
        panel.grid.minor.x=element_blank(), 
        panel.grid.minor.y=element_blank(),
        panel.grid.major.x=element_blank(), 
        panel.grid.major.y=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.ticks.x=element_blank(),
        axis.ticks.y=element_blank()) +
  annotate('rect', fill=panelcols[5], xmin=0.07, xmax=1, ymin=0.12, ymax=1, alpha=.5) +
  annotate('text', x=.55, y=.85, label='Strongman') +
  annotate('text', x=.55, y=.6, label='N = 615 (12%)', size=2.5) +
  annotate('text', x=.55, y=.3, label='MID country-years:\n126 (20%)', size=2.5)

## Runs models with defined scenarios

# Negative binomial by regime type, lagged MID = 0
predValueLag0 = exp(data.matrix(scenRegimesLag0) %*% t(negbdraws))
tomelt = t(predValueLag0)
colnames(tomelt) = c('Democracy', IVs[1:4])
l0nb = melt(tomelt, id.var=c())[,2:3]
colnames(l0nb) = c('regime', 'MIDs')
l0nb$regime = factor(l0nb$regime, levels=rev(c('Democracy', IVs[1:4])))

# Negative binomial by regime type, lagged MID = 1
predValueLag1 = exp(data.matrix(scenRegimesLag1) %*% t(negbdraws))
tomelt = t(predValueLag1)
colnames(tomelt) = c('Democracy', IVs[1:4])
l1nb = melt(tomelt, id.var=c())[,2:3]
colnames(l1nb) = c('regime', 'MIDs')
l1nb$regime = factor(l1nb$regime, levels=rev(c('Democracy', IVs[1:4])))

# Zero-inflated negative binomial by regime type, lagged MID = 0
predZerosLag0 = 1 / (1 + exp(-data.matrix(scenRegimesLag0) %*% t(zinbdraws[,11:20])))
predCountLag0 = exp(data.matrix(scenRegimesLag0) %*% t(zinbdraws[,1:10]))
predMMIDsLag0 = (1 - predZerosLag0) * predCountLag0
tomelt = t(predMMIDsLag0)
colnames(tomelt) = c('Democracy', IVs[1:4])
l0zi = melt(tomelt, id.var=c())[,2:3]
colnames(l0zi) = c('regime', 'MIDs')
l0zi$regime = factor(l0zi$regime, levels=rev(c('Democracy', IVs[1:4])))

# Zero-inflated negative binomial by regime type, lagged MID = 1
predZerosLag1 = 1 / (1 + exp(-data.matrix(scenRegimesLag1) %*% t(zinbdraws[,11:20])))
predCountLag1 = exp(data.matrix(scenRegimesLag1) %*% t(zinbdraws[,1:10]))
predMMIDsLag1 = (1 - predZerosLag1) * predCountLag1
tomelt = t(predMMIDsLag1)
colnames(tomelt) = c('Democracy', IVs[1:4])
l1zi = melt(tomelt, id.var=c())[,2:3]
colnames(l1zi) = c('regime', 'MIDs')
l1zi$regime = factor(l1zi$regime, levels=rev(c('Democracy', IVs[1:4])))

# Creates scenario-based density plots
densLag0 = list()
i = 1 # i must be manually incremented because looping produces blank plots
tmp = ggplot(data=subset(l0nb, regime %in% combn(c('Democracy', IVs[1:4]), 2)[,i]), 
             aes(x=MIDs, color=regime, fill=regime), position='dodge') + 
  geom_density(alpha=.5) +
  geom_density(data=subset(l0zi, regime %in% combn(c('Democracy', IVs[1:4]), 2)[,i]),
               alpha=.5, aes(y=-..density..)) +
  scale_color_manual(values=colcom[(2*i-1):(2*i)]) +
  scale_fill_manual(values=colcom[(2*i-1):(2*i)]) +
  xlab('') +
  ylab('') +
  guides(color=F, fill=F) +
  theme(panel.grid.minor.y=element_blank(), 
        panel.grid.major.y=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank())
tmp
densLag0[[i]] = tmp
i = i + 1
save(file='densLag0.RData', densLag0)

densLag1 = list()
i = 1 # i must be manually incremented because looping produces blank plots
tmp = ggplot(data=subset(l1nb, regime %in% combn(c('Democracy', IVs[1:4]), 2)[,i]), 
             aes(x=MIDs, color=regime, fill=regime), position='dodge') + 
  geom_density(alpha=.5) +
  geom_density(data=subset(l1zi, regime %in% combn(c('Democracy', IVs[1:4]), 2)[,i]),
               alpha=.5, aes(y=-..density..)) +
  scale_color_manual(values=colcom[(2*i-1):(2*i)]) +
  scale_fill_manual(values=colcom[(2*i-1):(2*i)]) +
  xlab('') +
  ylab('') +
  guides(color=F, fill=F) +
  theme(panel.grid.minor.y=element_blank(), 
        panel.grid.major.y=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank())
tmp
densLag1[[i]] = tmp
i = i + 1
save(file='densLag1.RData', densLag1)

# Plots all 20 scenario combinations with labels in diagonal
dens.L0.L1 = arrangeGrob(panel.dem, densLag0[[1]], densLag0[[2]], densLag0[[3]], densLag0[[4]],
                         densLag1[[1]], panel.mac, densLag0[[5]], densLag0[[6]], densLag0[[7]],
                         densLag1[[2]], densLag1[[5]], panel.jun, densLag0[[8]], densLag0[[9]],
                         densLag1[[3]], densLag1[[6]], densLag1[[8]], panel.bos, densLag0[[10]],
                         densLag1[[4]], densLag1[[7]], densLag1[[9]], densLag1[[10]], panel.str,
                         nrow=5)
save(file='densL0L1.RData', dens.L0.L1)
ggsave('densLag0vLag1v2.pdf', dens.L0.L1, width=9, height=9, units='in')

save(file='panel4_simulated_scenarios.RData', scenRegimesLag0, scenRegimesLag1, panelcols, 
     colcom, panel.dem, panel.mac, panel.jun, panel.bos, panel.str, l0nb, l0zi, l1nb, l1zi,
     densLag0, densLag1, dens.L0.L1)